The purpose of this Jupyter Notebook is to guides you through the process of Digital Elevation Model processing using open source tools. The case study focus on FKKN airspace, but the same process can be used for any other airspace.
Digital Elevation Model (DEM) is simply a mathematical representation of the continuous surface of the ground, based on a (large) number of points defined in terms of X, Y and Z co-ordinates. Typically, DEM data can be represented as a raster which is most easily expressed as being a 2D array with each individual cell having an elevation associated with it.
Raster or “gridded” data are stored as a grid of values which are rendered on a map as pixels. Each pixel value represents an area on the Earth’s surface. A raster file is a composed of regular grid of cells, all of which are the same size. You’ve looked at and used rasters before if you’ve looked at photographs or imagery in a tool like Google Earth. However, the raster files that you will work with are different from photographs in that they are spatially referenced. Each pixel represents an area of land on the ground. That area is defined by the spatial resolution of the raster.
Terrain and obstacles data have many applications in aviation insdustry:
The structure of FKKN airspace is describe in AIS ASECNA. We have three main features:
# Set path to FKKN TMA shapefile
fkkn_tma_path = os.path.join("fkkn_tma", "Fkkn_tma.shp")
fkkn_tma = gpd.read_file(fkkn_tma_path)
# Set path to FKKN CTR shapefile
fkkn_ctr_path = os.path.join("fkkn_ctr", "Fkkn_ctr.shp")
fkkn_ctr = gpd.read_file(fkkn_ctr_path)
# Set path to FKP3 shapefile
fkkn_kp3_path = os.path.join("fkkn_fkp3", "Fkkn_fkp3.shp")
fkkn_fkp3 = gpd.read_file(fkkn_kp3_path)
fkkn_tma["Name"]
0 TMA FKKN Name: Name, dtype: object
fig, ax = plt.subplots(figsize=(20, 20))
# Plot FKKN TMA
fkkn_tma.plot( legend= True,
figsize=(10, 6),
cmap="Set2",
edgecolor="black",
facecolor= "none",
linewidth=3,
#linestyle='dotted',
ax=ax)
# Plot FKKN CTR
fkkn_ctr.plot(
edgecolor="blue",
facecolor="none",
linewidth=3,
legend=True,
linestyle='dotted',
ax=ax)
# Plot FKP3
fkkn_fkp3.plot(
edgecolor="red",
facecolor= None,
color = "white",
hatch="\\",
linewidth=3,
legend = True,
ax=ax)
"""
ax.set(xlabel="LONGITUDE ( Decimal Degrees)",
ylabel="LATITUDE (Decimal Degrees)",
title="FKKN AIRSPACE")
"""
ax.set_title("FKKN AIRSPACE", color = "black",fontsize=20)
#ax.set_xlabel("Longitude ( Decimal Degrees)",color = "blue",style='italic',fontsize=16, fontweight='bold')
ax.set_xlabel("LONGITUDE ( Decimal Degrees)",color = "black",fontsize=18)
ax.set_ylabel("LATITUDE (Decimal Degrees)",color = "black",fontsize=18)
ax.set_axisbelow(True)
ax.yaxis.grid(color='gray', linestyle='dashed')
ax.xaxis.grid(color='gray', linestyle='dashed')
props = dict(boxstyle='round', edgecolor='black', facecolor = "white", alpha=0.2)
ax.text(0.45, 0.90, 'FKKN TMA' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='black', fontweight='bold', )
ax.text(0.40, 0.45, 'CTR FKKN' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='blue', fontweight='bold', )
ax.text(0.30, 0.65, 'FKP3' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='red', fontweight='bold', )
plt.savefig('01_fkkn_tma.png')
plt.show()
#ctx.add_basemap(ax, url = ctx.providers.OpenStreetMap.Mapnik)
We use OpenStreetMap.Mapnik base map
fig, ax = plt.subplots(figsize=(20, 20))
# Plot FKKN TMA
fkkn_tma.plot( legend= True,
figsize=(10, 6),
cmap="Set2",
edgecolor="black",
facecolor= "none",
linewidth=3,
#linestyle='dotted',
ax=ax)
# Plot FKKN CTR
fkkn_ctr.plot(
edgecolor="blue",
facecolor="none",
linewidth=3,
legend=True,
linestyle='dotted',
ax=ax)
# Plot FKP3
fkkn_fkp3.plot(
edgecolor="red",
facecolor= "None",
color = "None",
hatch="\\",
linewidth=3,
legend = True,
alpha=0.5,
ax=ax)
"""
ax.set(xlabel="Longitude ( Decimal Degrees)",
ylabel="Latitude (Decimal Degrees)",
title="FKKN AIRSPACE WITH MAPNIK BASEMAP")
"""
ax.set_title("FKKN AIRSPACE WITH MAPNIK BASEMAP", color = "black",fontsize=20)
ax.set_xlabel("LONGITUDE ( Decimal Degrees)",color = "black",fontsize=18)
ax.set_ylabel("LATITUDE (Decimal Degrees)",color = "black",fontsize=18)
ax.set_axisbelow(True)
ax.yaxis.grid(color='gray', linestyle='dashed')
ax.xaxis.grid(color='gray', linestyle='dashed')
ctx.add_basemap(ax, crs=fkkn_tma.crs.to_string(), source=ctx.providers.OpenStreetMap.Mapnik)
#ctx.add_basemap(ax, crs=fkkn_tma.crs.to_string(), source=ctx.providers.Esri.WorldPhysical)
props = dict(boxstyle='round', edgecolor='black', facecolor = "white", alpha=0.2)
ax.text(0.45, 0.90, 'FKKN TMA' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='black', fontweight='bold', )
ax.text(0.40, 0.45, 'CTR FKKN' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='blue', fontweight='bold', )
ax.text(0.30, 0.65, 'FKP3' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='red', fontweight='bold', )
plt.savefig('02_fkkn_tma.png')
plt.show()
#ctx.add_basemap(ax, crs=fkkn_tma.crs.to_string(), source=ctx.providers.Esri.WorldPhysical)
The following statement creates an buffer airspace of 10 NM around the the FKKN TMA.
We choose 1ONM buffuring because, we will use that portion of airspace later for Mountainous area calculation inside TMA FKKN in accordance with ICAO provisions in DOC DOC 8168 VOL II, PANS OBS
from qgis.analysis import QgsNativeAlgorithms
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
import processing
#The path to TMA FKKN file as input
intput_path = os.path.join("fkkn_tma", "Fkkn_tma.shp")
#The path to TMA FKKN + 1ONM BUFFER file as output
output_path = os.path.join("fkkn_tma", "fkkn_tma_buffer_10NM.shp")
result = processing.run("native:buffer",
{
'INPUT':intput_path,
'DISTANCE':.1667,
#'SEGMENTS':5,
'END_CAP_STYLE':0,
'JOIN_STYLE':0,
#'MITER_LIMIT':2,
'DISSOLVE':False,
'OUTPUT':output_path,
}
)
buffer = result['OUTPUT']
Let see the buffer we just create
# load the buffer shapefile for visualisation
output_path = os.path.join("fkkn_tma", "fkkn_tma_buffer_10NM.shp")
fkkn_tma_buffer = gpd.read_file(output_path)
fig, ax = plt.subplots(figsize=(20, 20))
# Plot the data and add a legend
fkkn_tma_buffer.plot(
legend=True,
figsize=(10, 6),
#cmap="Set2",
edgecolor="green",
facecolor="none",
linewidth=3,
ax=ax)
# Plot FKKN TMA
fkkn_tma.plot( legend= True,
figsize=(10, 6),
cmap="Set2",
edgecolor="black",
facecolor= "none",
linewidth=3,
#linestyle='dotted',
ax=ax)
# Plot FKKN CTR
fkkn_ctr.plot(
edgecolor="blue",
facecolor="none",
linewidth=3,
legend=True,
linestyle='dotted',
ax=ax)
# Plot FKP3
fkkn_fkp3.plot(
edgecolor="red",
facecolor= None,
color = "white",
hatch="\\",
linewidth=3,
legend = True,
#alpha=0.5,
ax=ax)
"""
ax.set(xlabel="Longitude ( Decimal Degrees)",
ylabel="Latitude (Decimal Degrees)",
title="FKKN AIRSPACE + 10NM BUFFER")
"""
ax.set_title("FKKN AIRSPACE + 10NM BUFFER", color = "black",fontsize=20)
ax.set_xlabel("LONGITUDE ( Decimal Degrees)",color = "black",fontsize=18)
ax.set_ylabel("LATITUDE (Decimal Degrees)",color = "black",fontsize=18)
ax.set_axisbelow(True)
ax.yaxis.grid(color='gray', linestyle='dashed')
ax.xaxis.grid(color='gray', linestyle='dashed')
#ctx.add_basemap(ax, crs=fkkn_tma.crs.to_string(), source=ctx.providers.OpenStreetMap.Mapnik)
#ctx.add_basemap(ax, crs=fkkn_tma.crs.to_string(), source=ctx.providers.Esri.WorldPhysical)
props = dict(boxstyle='round', edgecolor='black', facecolor = "white", alpha=0.2)
ax.text(0.40, 0.93, '10NM BUFFER AROUND FKKN TMA' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='green', fontweight='bold', )
ax.text(0.45, 0.80, 'FKKN TMA' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='black', fontweight='bold', )
ax.text(0.40, 0.45, 'CTR FKKN' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='blue', fontweight='bold', )
ax.text(0.37, 0.63, 'FKP3' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='red', fontweight='bold', )
plt.savefig('03_fkkn_tma.png')
plt.show()
The DEM we will process in this section is a lager data set covering the entire Cameroon.
It is provided by Cameroon Civil Aviation Authority ATS Bureau.
The same DEM is used to design Instrument Approach Procedure of Ngaoundéré Aerodrome.
liste = !ls cmr_dem_folder/Data
print(len(liste))
268
print(liste)
['e008n03.BIL', 'e008n03.hdr', 'e008n03.prj', 'e008n03.stx', 'e008n04.BIL', 'e008n04.hdr', 'e008n04.prj', 'e008n04.stx', 'e008n05.BIL', 'e008n05.hdr', 'e008n05.prj', 'e008n05.stx', 'e009n02.BIL', 'e009n02.hdr', 'e009n02.prj', 'e009n02.stx', 'e009n03.BIL', 'e009n03.hdr', 'e009n03.prj', 'e009n03.stx', 'e009n04.BIL', 'e009n04.hdr', 'e009n04.prj', 'e009n04.stx', 'e009n05.BIL', 'e009n05.hdr', 'e009n05.prj', 'e009n05.stx', 'e009n06.BIL', 'e009n06.hdr', 'e009n06.prj', 'e009n06.stx', 'e010n02.BIL', 'e010n02.hdr', 'e010n02.prj', 'e010n02.stx', 'e010n03.BIL', 'e010n03.hdr', 'e010n03.prj', 'e010n03.stx', 'e010n04.BIL', 'e010n04.hdr', 'e010n04.prj', 'e010n04.stx', 'e010n05.BIL', 'e010n05.hdr', 'e010n05.prj', 'e010n05.stx', 'e010n06.BIL', 'e010n06.hdr', 'e010n06.prj', 'e010n06.stx', 'e010n07.BIL', 'e010n07.hdr', 'e010n07.prj', 'e010n07.stx', 'e011n02.BIL', 'e011n02.hdr', 'e011n02.prj', 'e011n02.stx', 'e011n03.BIL', 'e011n03.hdr', 'e011n03.prj', 'e011n03.stx', 'e011n04.BIL', 'e011n04.hdr', 'e011n04.prj', 'e011n04.stx', 'e011n05.BIL', 'e011n05.hdr', 'e011n05.prj', 'e011n05.stx', 'e011n06.BIL', 'e011n06.hdr', 'e011n06.prj', 'e011n06.stx', 'e011n07.BIL', 'e011n07.hdr', 'e011n07.prj', 'e011n07.stx', 'e011n08.BIL', 'e011n08.hdr', 'e011n08.prj', 'e011n08.stx', 'e012n02.BIL', 'e012n02.hdr', 'e012n02.prj', 'e012n02.stx', 'e012n03.BIL', 'e012n03.hdr', 'e012n03.prj', 'e012n03.stx', 'e012n04.BIL', 'e012n04.hdr', 'e012n04.prj', 'e012n04.stx', 'e012n05.BIL', 'e012n05.hdr', 'e012n05.prj', 'e012n05.stx', 'e012n06.BIL', 'e012n06.hdr', 'e012n06.prj', 'e012n06.stx', 'e012n07.BIL', 'e012n07.hdr', 'e012n07.prj', 'e012n07.stx', 'e012n08.BIL', 'e012n08.hdr', 'e012n08.prj', 'e012n08.stx', 'e012n09.BIL', 'e012n09.hdr', 'e012n09.prj', 'e012n09.stx', 'e013n02.BIL', 'e013n02.hdr', 'e013n02.prj', 'e013n02.stx', 'e013n03.BIL', 'e013n03.hdr', 'e013n03.prj', 'e013n03.stx', 'e013n04.BIL', 'e013n04.hdr', 'e013n04.prj', 'e013n04.stx', 'e013n05.BIL', 'e013n05.hdr', 'e013n05.prj', 'e013n05.stx', 'e013n06.BIL', 'e013n06.hdr', 'e013n06.prj', 'e013n06.stx', 'e013n07.BIL', 'e013n07.hdr', 'e013n07.prj', 'e013n07.stx', 'e013n08.BIL', 'e013n08.hdr', 'e013n08.prj', 'e013n08.stx', 'e013n09.BIL', 'e013n09.hdr', 'e013n09.prj', 'e013n09.stx', 'e013n10.BIL', 'e013n10.hdr', 'e013n10.prj', 'e013n10.stx', 'e013n11.BIL', 'e013n11.hdr', 'e013n11.prj', 'e013n11.stx', 'e013n12.BIL', 'e013n12.hdr', 'e013n12.prj', 'e013n12.stx', 'e014n01.BIL', 'e014n01.hdr', 'e014n01.prj', 'e014n01.stx', 'e014n02.BIL', 'e014n02.hdr', 'e014n02.prj', 'e014n02.stx', 'e014n03.BIL', 'e014n03.hdr', 'e014n03.prj', 'e014n03.stx', 'e014n04.BIL', 'e014n04.hdr', 'e014n04.prj', 'e014n04.stx', 'e014n05.BIL', 'e014n05.hdr', 'e014n05.prj', 'e014n05.stx', 'e014n06.BIL', 'e014n06.hdr', 'e014n06.prj', 'e014n06.stx', 'e014n07.BIL', 'e014n07.hdr', 'e014n07.prj', 'e014n07.stx', 'e014n08.BIL', 'e014n08.hdr', 'e014n08.prj', 'e014n08.stx', 'e014n09.BIL', 'e014n09.hdr', 'e014n09.prj', 'e014n09.stx', 'e014n10.BIL', 'e014n10.hdr', 'e014n10.prj', 'e014n10.stx', 'e014n11.BIL', 'e014n11.hdr', 'e014n11.prj', 'e014n11.stx', 'e014n12.BIL', 'e014n12.hdr', 'e014n12.prj', 'e014n12.stx', 'e014n13.BIL', 'e014n13.hdr', 'e014n13.prj', 'e014n13.stx', 'e015n01.BIL', 'e015n01.hdr', 'e015n01.prj', 'e015n01.stx', 'e015n02.BIL', 'e015n02.hdr', 'e015n02.prj', 'e015n02.stx', 'e015n03.BIL', 'e015n03.hdr', 'e015n03.prj', 'e015n03.stx', 'e015n04.BIL', 'e015n04.hdr', 'e015n04.prj', 'e015n04.stx', 'e015n06.BIL', 'e015n06.hdr', 'e015n06.prj', 'e015n06.stx', 'e015n07.BIL', 'e015n07.hdr', 'e015n07.prj', 'e015n07.stx', 'e015n08.BIL', 'e015n08.hdr', 'e015n08.prj', 'e015n08.stx', 'e015n09.BIL', 'e015n09.hdr', 'e015n09.prj', 'e015n09.stx', 'e015n10.BIL', 'e015n10.hdr', 'e015n10.prj', 'e015n10.stx', 'e015n11.BIL', 'e015n11.hdr', 'e015n11.prj', 'e015n11.stx', 'e015n12.BIL', 'e015n12.hdr', 'e015n12.prj', 'e015n12.stx', 'e016n01.BIL', 'e016n01.hdr', 'e016n01.prj', 'e016n01.stx', 'e016n02.BIL', 'e016n02.hdr', 'e016n02.prj', 'e016n02.stx', 'e016n03.BIL', 'e016n03.hdr', 'e016n03.prj', 'e016n03.stx']
We note that, the folder containing the entire DEM have 268 files representing 67 blocks of DEM.
The data format for each block is a Band Interleaved by Line: .BIL.
The meta data are depicted in 3 differents files with extension:
.hdr : header informmation.proj : projection information.stx: statistics information (band number, minimum value, maximum value, mean value, and standard deviation)Open a rondom block for data exploration
dem_path = os.path.join("cmr_dem_folder","Data", "e012n08.BIL")
dem = rio.open(dem_path)
DEM coordinate reference sustem
print(dem.crs)
EPSG:4326
DEM coordinate reference sustem as wel known text (wkt)
print(dem.crs.to_wkt())
GEOGCS["Geographic Coordinate System",DATUM["WGS84",SPHEROID["WGS84",6378137,298.257223560493]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]
DEM Resolution
print(dem.res)
(0.000833333333, 0.000833333333)
Number of Bands
print(dem.count)
1
The width
print(dem.width)
1201
The height
print(dem.height)
1201
The Driver (data format)
print(dem.driver)
EHdr
No data values for all channels
print(dem.nodatavals)
(-9999.0,)
All Metadata for the whole raster dataset
print(dem.meta)
{'driver': 'EHdr', 'dtype': 'int16', 'nodata': -9999.0, 'width': 1201, 'height': 1201, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.000833333333, 0.0, 11.9995833333335,
0.0, -0.000833333333, 9.0004166666665)}
We saw previously, that our DEM is split into 67 blocks of .BIL file. In this section , we will merge all small blocks of DEM into a single larger Geotiff file cameroon.tif.
The challenge here is to do the operation without loosing any data or meta data information
! ls *.BIL > filelist.txt
#merge all small DEM block in a virtual raster layer : cameroon.vrt
! gdalbuildvrt -input_file_list filelist.txt cameroon.vrt
0...10...20...30...40...50...60...70...80...90...100 - done.
#Translate the virtual raster layer (cameroon.vrt) to a Geotiff raster format (cameroon.tif)
! gdal_translate -of GTiff cameroon.vrt cameroon.tif
Input file size is 10801, 15601 0...10...20...30...40...50...60...70...80...90...100 - done.
#Load the raster file we just created
dem_cmr = rio.open('cameroon.tif')
#print the cordinate reference system
print(dem_cmr.crs)
EPSG:4326
#print the cordinate reference system as WKT
print(dem_cmr.crs.to_wkt())
GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.25722356049,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
#print the resolution
print(dem_cmr.res)
(0.0008333333330000007, 0.0008333333330000007)
#print the numbers of band
print(dem_cmr.count)
1
#Nodata values
print(dem_cmr.nodatavals)
(-9999.0,)
#print meta data information
print(dem_cmr.meta)
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -9999.0, 'width': 10801, 'height': 15601, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.0008333333330000007, 0.0, 7.9995833333335,
0.0, -0.0008333333330000007, 14.0004166666665)}
We note that all these information is the same that we got before processing. We did not loss information during the merging processs
As we see previously, the DEM covers the entire cameroon, it means more than 474 444 square kilometers. We are highly motivated to reduce the size of an image to only our area of interest (AOI). The best way to accomplish this reduction is to clip an image to a boundary that defines our study area.
Let see how look the DEM that we want to clip.
cmr_dem_path = os.path.join("cmr_dem_folder","Data", "cameroon.tif")
# open the lidar chm
with rio.open(cmr_dem_path ) as src:
cmr_dem_im = src.read(masked=True)[0]
extent = rio.plot.plotting_extent(src)
cmr_dem_profile = src.profile
fig, ax = plt.subplots(figsize=(20, 20))
ep.plot_bands(cmr_dem_im ,
#cmap='terrain',
#cmap = 'magma',
cmap='viridis_r',
#cmap="jet",
extent=extent,
title="DEM COVERING CAMEROON BOUNDARY",
cbar = True,
ax = ax,
alpha = 0.5)
plt.savefig('04_fkkn_tma.png')
plt.show()
#plt.savefig('04_fkkn_tma.png')
cmr_dem_path = os.path.join("cmr_dem_folder","Data", "cameroon.tif")
cmr_dem = rio.open(cmr_dem_path)
fig, ax = plt.subplots(figsize=(20, 20))
"""
show(cmr_dem,
cmap='viridis_r',
#cmap='Greys',
alpha=0.5,
#cmap='viridis',
#title="Cameroon Digital Elevation Model (DEM) \n Inside FKKN TMA airspace + 10NM Buffer",
ax=ax)
"""
ep.plot_bands(cmr_dem_im ,
#cmap='terrain',
#cmap = 'magma',
cmap='viridis_r',
#cmap="jet",
extent=extent,
title="Cameroon Digital Elevation Model(DEM)",
cbar = True,
ax = ax)
fkkn_tma.plot(
legend=True,
#cmap="Set2",
edgecolor="magenta",
facecolor="None",
linewidth=3,
ax=ax)
"""
ax.set_axis_off()
plt.savefig('05_fkkn_tma.png')
plt.show()
"""
ax.set_title("DEM COVERING CAMEROON BOUNDARY \n INSIDE \n FKKN TMA + 10 NM BUFFER", color = "black",fontsize=20)
#ax.set_xlabel("LONGITUDE ( Decimal Degrees)",color = "black",fontsize=18)
#ax.set_ylabel("LATITUDE (Decimal Degrees)",color = "black",fontsize=18)
#ax.set_axisbelow(True)
#ax.yaxis.grid(color='gray', linestyle='dashed')
#ax.xaxis.grid(color='gray', linestyle='dashed')
plt.savefig('05_fkkn_tma.png')
plt.show()
cmr_dem.close()
We extract our AOI in color : magenta on the image above
#The path to our DEM
#input_dem = os.path.join('cmr_dem_folder', 'data', 'cameroon.tif')
input_dem= os.path.join("cmr_dem_folder","Data", "cameroon.tif")
#The path to TMA FKKN + 1ONM BUFFER fto be used as mask for clipping
mask_path = os.path.join("fkkn_tma", "fkkn_tma_buffer_10NM.shp")
#The path to TMA FKKN + 1ONM BUFFER DEM as output
output_dem = os.path.join('cmr_dem_folder', 'Data', 'dem_fkkn_tma_buffer_10NM.tif')
result = processing.run( "gdal:cliprasterbymasklayer",
{
'INPUT':input_dem,
'MASK': mask_path,
'OUTPUT':output_dem,
}
)
clipped_dem = result['OUTPUT']
dem_fkkn_tma_buffer_10NM_path = os.path.join('cmr_dem_folder', 'Data', 'dem_fkkn_tma_buffer_10NM.tif')
with rio.open(dem_fkkn_tma_buffer_10NM_path) as src:
dem_fkkn_tma_buffer_10NM = src.read(1, masked=True)
dem_ext = rio.plot.plotting_extent(src)
fig, ax = plt.subplots(figsize=(20, 20))
ep.plot_bands(dem_fkkn_tma_buffer_10NM,
cmap='jet',
#cmap='Greys',
extent=dem_ext,
#title="DEM MODEL COVERING FKKN TMA + 10NM BUFFER",
cbar=True,
alpha = 0.5,
ax = ax)
ax.set_title("DEM MODEL COVERING FKKN TMA + 10NM BUFFER", color = "black",fontsize=20)
#ax.set_xlabel("LONGITUDE ( Decimal Degrees)",color = "black",fontsize=18)
#ax.set_ylabel("LATITUDE (Decimal Degrees)",color = "black",fontsize=18)
#ax.set_axisbelow(True)
#ax.yaxis.grid(color='gray', linestyle='dashed')
#ax.xaxis.grid(color='gray', linestyle='dashed')
plt.savefig('06_fkkn_tma.png')
plt.show()
#plt.savefig('04_fkkn_tma.png')
dem_fkkn_tma_buffer_10NM = rio.open(dem_fkkn_tma_buffer_10NM_path)
fig, ax = plt.subplots(figsize=(20, 20))
ep.plot_bands(dem_fkkn_tma_buffer_10NM,
cmap='jet',
#cmap='Greys',
extent=dem_ext,
title="DEM MODEL COVERING FKKN TMA + 10NM BUFFER",
cbar=True,
alpha = 0.5,
ax = ax)
# Plot the data and add a legend
fkkn_tma_buffer.plot(
legend=True,
#figsize=(10, 6),
#cmap="Set2",
edgecolor="magenta",
facecolor="none",
linewidth=3,
ax=ax)
# Plot FKKN TMA
fkkn_tma.plot( legend= True,
#figsize=(10, 6),
cmap="Set2",
edgecolor="black",
facecolor= "none",
linewidth=3,
#linestyle='dotted',
ax=ax)
# Plot FKKN CTR
fkkn_ctr.plot(
edgecolor="blue",
facecolor="none",
linewidth=3,
legend=True,
linestyle='dotted',
ax=ax)
# Plot FKP3
fkkn_fkp3.plot(
edgecolor="red",
facecolor= "None",
#color = "white",
hatch="\\",
linewidth=3,
legend = True,
#alpha=0.5,
ax=ax)
"""
ax.set(xlabel="Longitude ( Decimal Degrees)",
ylabel="Latitude (Decimal Degrees)",
title="DEM COVERING FKKN AIRSPACE + 10 NM BUFFER")
"""
ax.set_title("DEM COVERING FKKN TMA + 10NM \n INSIDE \n FKKN AIRSPACE", color = "black",fontsize=20)
ax.set_xlabel("LONGITUDE ( Decimal Degrees)",color = "black",fontsize=18)
ax.set_ylabel("LATITUDE (Decimal Degrees)",color = "black",fontsize=18)
ax.set_axisbelow(True)
ax.yaxis.grid(color='gray', linestyle='dashed')
ax.xaxis.grid(color='gray', linestyle='dashed')
props = dict(boxstyle='round', edgecolor='black', facecolor = "None", alpha=0.2)
ax.text(0.30, 0.93, '10NM BUFFER AROUND FKKN TMA' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='red', fontweight='bold', )
ax.text(0.40, 0.80, 'FKKN TMA' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='red', fontweight='bold', )
ax.text(0.40, 0.45, 'CTR FKKN' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='red', fontweight='bold', )
ax.text(0.33, 0.63, 'FKP3' , transform=ax.transAxes, fontsize=16, verticalalignment='top',bbox=props, color='red', fontweight='bold', )
plt.savefig('07_fkkn_tma.png')
plt.show()
Now that we have extrated the DEM of our AOI, It is time to do some analysis and products derivations.
In this part, we will cover:
When we apply a ramp of color to elevation, it gives us a quick overview of how our elevation data look like. It is easy to know where we have high/low terrain
from PyQt5.QtGui import QColor
data_dir = os.path.join('cmr_dem_folder', 'Data')
raster_name = 'dem_fkkn_tma_buffer_10NM.tif'
raster_path = os.path.join(data_dir, raster_name)
lyr = QgsRasterLayer(raster_path, "DEM")
s = QgsRasterShader()
c = QgsColorRampShader()
i = []
i.append(c.ColorRampItem(300, QColor('#d7191c'), '1000 feets'))
i.append(c.ColorRampItem(600, QColor('#fdae61'), '2000 feets'))
i.append(c.ColorRampItem(900, QColor('#ffffbf'), '3000 feets'))
i.append(c.ColorRampItem(1200, QColor('#abdda4'), '4000 feet'))
i.append(c.ColorRampItem(1500, QColor('#2b83ba'), '5000 feets'))
c.setColorRampItemList(i)
s.setRasterShaderFunction(c)
ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), 1, s)
lyr.setRenderer(ps)
#QgsMapLayerRegistry.instance().addMapLayer(lyr)
lyr.isValid()
QgsProject.instance().addMapLayer(lyr)
from IPython.display import Image
Image("color_ramp.png")
Shaded relief maps color elevation in such a way that it looks as if the terrain is cast in a low angle light, which creates bright spots and shadows. This aesthetic styling creates an almost photographic illusion, which is easy to grasp so that we can understand the variation in the terrain. It is important to note that this style is truly an illusion as the light is often physically inaccurate in terms of the solar angle, and the elevation is usually exaggerated to increase contrast.
import os
import processing
data_dir = os.path.join('cmr_dem_folder', 'Data')
raster_name = 'dem_fkkn_tma_buffer_10NM.tif'
raster_path = os.path.join(data_dir, raster_name)
raster_out_name = "dem_fkkn_tma_buffer_hillshade.tif"
raster_out_path = os.path.join(data_dir, raster_out_name)
results = processing.runAndLoadResults("native:hillshade",
{'INPUT': raster_path,
'Z_FACTOR':3, # Elevation exageration: 1 mean true scale
'AZIMUTH':315,# Sun direction: between 0 and 360 degrees
'V_ANGLE':45, # Sun angle: between 0 and 90 degres
'OUTPUT': raster_out_path})
data_dir = os.path.join('cmr_dem_folder', 'Data')
raster_out_name = "dem_fkkn_tma_buffer_hillshade.tif"
raster_out_path = os.path.join(data_dir, raster_out_name)
hillshade_dem = rio.open(raster_out_path)
fig, ax = plt.subplots(figsize=(20, 10))
show(hillshade_dem ,
cmap='viridis_r',
#cmap='Greys',
alpha=0.5,
#cmap='viridis',
title="FKKN TMA DEM HILLS SHADE",
ax=ax)
ax.set_axis_off()
A contour is an isoline along the same elevation in a dataset. Contours are usually stepped at intervals to create an intuitive way to represent elevation data, both visually and numerically, using a resource-efficient vector dataset. Now, let's look at another way to visualize the elevation better using contours.
import os
import processing
data_dir = os.path.join('cmr_dem_folder', 'Data')
raster_name = 'dem_fkkn_tma_buffer_10NM.tif'
raster_path = os.path.join(data_dir, raster_name)
data_dir2 = os.path.join('cmr_dem_folder', 'Data','fkkn_buffer_contour')
#
shapefile_name = "dem_fkkn_buffer_contours.shp"
shapefile_path = os.path.join(data_dir2, shapefile_name)
#
rasterLyr = QgsRasterLayer(raster_path, "DEM")
rasterLyr.isValid()
QgsProject.instance().addMapLayer(rasterLyr)
processing.runAndLoadResults("gdal:contour",
{'INPUT':raster_path,
'BAND':1,
'INTERVAL':300,# adjust interval as required
'FIELD_NAME':'ELEV',
'OUTPUT': shapefile_path})
data_dir2 = os.path.join('cmr_dem_folder', 'Data','fkkn_buffer_contour')
#
shapefile_name = "dem_fkkn_buffer_contours.shp"
shapefile_path = os.path.join(data_dir2, shapefile_name)
# load the buffer shapefile for visualisation
fkkn_tma_buffer_contour = gpd.read_file(shapefile_path )
fig, ax = plt.subplots(figsize=(20, 15))
# Plot FKKN TMA
fkkn_tma_buffer_contour.plot(legend=True,
#figsize=(10, 6),
cmap="Set2",
edgecolor="magenta",
facecolor="none",
linewidth=0.2,
ax=ax)
ax.set_axis_off()
ax.set(#xlabel="Longitude ( Decimal Degrees)",
#ylabel="Latitude (Decimal Degrees)",
title="FKKN DEM ELEVATION CONTOUR LINE ")
[Text(0.5, 1.0, 'FKKN DEM ELEVATION CONTOUR LINE ')]
The image above show the elevation contour for the entire AOI, to better understand what going on, let zoom in into a small area, to see how it look like
Image("dem_contour.png")
As mentionned at the begining of this work, DEM is used by Instrument Flight Procédure Designer (IFPD) for many tasks including (MOC calculation, Mountainous Area calculation, MSA, MDA/H or DA/DH, Contingency procedures, Drift-Down Procedures, Emergency En-route Landing, etc.).
DEM model is also used as background for many aeronautical charts like (Instruments Approach Chart(IAC), Standard Arrivals (STARs), Standard Instruments Departure (SID), etc.)
The final result of our analysis, is a product that can be used to design an Instrument procecure of Ngaoundéré Airport. It is important to mention that, the raw/original data that we use, is the same that CCAA uses for the ongoing process of Ngaoundere procedure design.
For example, aeronautical charts provide by JEPPESEN are georeferenced(each pixel on the chart match a physical location on the ground). When use in the onboard application like foreflight, the DEM model on the background of chart can significantly inscrease pilot awareness with respect to terrain clearance in a real time.
We appliy color ramp to elevation contour line to get a DEM that can be use as a background map for aeronautical chart
from PyQt5.QtGui import QColor
data_dir = os.path.join('cmr_dem_folder', 'Data')
raster_name = 'dem_fkkn_tma_buffer_10NM.tif'
raster_path = os.path.join(data_dir, raster_name)
lyr = QgsRasterLayer(raster_path, "DEM")
s = QgsRasterShader()
c = QgsColorRampShader()
i = []
i.append(c.ColorRampItem(300, QColor('#d7191c'), '1000 feets'))
i.append(c.ColorRampItem(600, QColor('#fdae61'), '2000 feets'))
i.append(c.ColorRampItem(900, QColor('#ffffbf'), '3000 feets'))
i.append(c.ColorRampItem(1200, QColor('#abdda4'), '4000 feet'))
i.append(c.ColorRampItem(1500, QColor('#2b83ba'), '5000 feets'))
c.setColorRampItemList(i)
s.setRasterShaderFunction(c)
ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), 1, s)
lyr.setRenderer(ps)
Image("final_product-dem_ynd_final.png")
Let zoom in to see how it look like
Image("contour_color_ramp.png")
The picture below show how elevation contours line and color ramp are use as background for aeronautical chart
Example 1: FKKL STARs RWY31
from IPython.display import Image
Image("fkkl_start_rwy13.png")
Example 2: FKKL IAC RNAV(GNSS)-RWY 31
from IPython.display import Image
Image("IAC_RWY31.png")
In this Jupyter notebbok, from the description of FKKN TMA in the AIS ASECNA, we build a spatial representation of this airspace in shapefile format. Then we use that shapefile to create a buffer of 10NM around . From a lager Digital Elevation Model(DEM) covering the entire Cameroon, we use that buffer as a mask
to extract the portion covering our Area of Interest(AOI).
After processing the DEM of our AOI, we get an Elevation Contours line and Color Ramp which can be used as background for aeronatical charts (IAC, STARs, SID, etc.).
We also apply Reproducible science technique to this work,so anyone can understand and replicate the steps of the analysis, applied to the same or even new data. We can use the next step to any other airpace and get the result.
The next release will cover Mountainous Area calculation for Instrument Flight Procedure Design (IFPD).
The Mountainous area. An area of changing terrain profile where the changes of terrain elevation exceed 900 m (3 000 ft) within a distance of 18.5 km (10.0 NM). DOC 8168 VOL II.
The work will consist of finding the points (lat/long) in WGS84 coordinate reference system where the difference between the highest elevation and the lower elevation within 10NM is 900m (3000ft) or above.
These points are important for MSA and MOC calculation, as we refer to Annexe 2,Chapter 5, alinea 5.1.2, Minimum Flight Altitude ( Except when necessary for take-off or landing, or except when specifically authorized by the appropriate authority) is 2000 ft over high terrain or in mountainous area and 1000 ft elsewhere, above the highest obstacle located within 8 km of the estimated position of the aircraft.
The intended output will be :